{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.2.1 RVM for regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "目标变量t的条件概率分布为:\n",
    "$$p(t|x,w\\beta)=\\mathcal N(t|y(x),\\beta^{-1})$$\n",
    "均值由线性模型给出:$$y(x)=w^T\\phi(x)$$\n",
    "写成SVM相类似的形式:\n",
    "$$y(x)=\\sum_{n=1}^Nw_nk(x,x_n)+b$$\n",
    "参数的数量为$M = N + 1$,$y(x)$与SVM的预测模型具有相同的形式,与SVM的情形相反，没有正定核的限制，基函数也没有被训练数据点的数量或位置所限制。\n",
    "\n",
    "似然函数为:\n",
    "$$p(t|X,w,\\beta)=\\prod^N_{n=1}p(t_n|x_n,w,\\beta)$$\n",
    "引入w的先验分布，为每一个$w_i$,都引入一个超参数$\\alpha_i$(ard先验):\n",
    "$$p(w|\\alpha)=\\prod_{i=1}^N\\mathcal N(w_i|0,\\alpha_i^{-1})$$\n",
    "当我们关于这些超参数\n",
    "最⼤化模型证据时，⼤部分都趋于⽆穷，对应的权参数的后验概率分布集中在零附近,⽣成了⼀个稀疏的模型。\n",
    "w的后验分布为:$$p(w|t,X,\\alpha,\\beta)=\\mathcal N(w|m,\\Sigma)$$\n",
    "均值为:$$m=\\beta\\Sigma\\Phi^Tt$$\n",
    "方差为:$$\\Sigma=(A+\\beta\\Phi^T\\Phi)^{-1}$$\n",
    "\n",
    "$\\Phi$是NxM的设计矩阵，\n",
    "元素为$\\Phi_{ni}=\\phi_i(x_n)(i=1,...,N),\\Phi_{nM}=1(n=1,...,N),A=diag(\\alpha_{i})$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "边缘似然函数为:\n",
    "$$p(t|X,\\alpha,\\beta)=\\int p(t|X,w,\\beta)p(w|\\alpha)dw$$\n",
    "对数边缘似然函数为:\n",
    "$$\\begin{align}lnp(t|X,\\alpha,\\beta)&=ln\\mathcal N(t|0,C)\\\\\n",
    "&=-\\frac{1}{2}\\{Nln(2\\pi)+ln|c|+t^TC^{-1}t\\}\\end{align}$$\n",
    "其中C为NXN矩阵:\n",
    "$$C=\\beta^{-1}I+\\Phi A^{-1}\\Phi^T$$\n",
    "\n",
    "关于超参数$\\alpha$,$\\beta$最大化边缘似然函数得到$\\alpha^*,\\beta^*$\n",
    "现在要关于$\\alpha$和$\\beta$最大化对数边缘似然函数，采用类似线性模型中的证据近似方法\n",
    "* $$\\alpha_i^{new}=\\frac{\\gamma_i}{m_i^2}$$\n",
    "* $$(\\beta^{new})^{-1}=\\frac{\\lVert t-\\Phi m\\lVert^2}{N-\\sum_i\\gamma_i}$$\n",
    "\n",
    "其中$m_i$是后验均值$m$的第i个分量，$\\gamma_i=1-\\alpha_i\\sum_{ii}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "预测分布为:\n",
    "$$\\begin{align}p(t|x,X,t,\\alpha^*,\\beta^*)&=\\int p(t|x,w,\\beta*)p(w|X,t,\\alpha^*,\\beta^*)dw\\\\&=\\mathcal N(t|m^T\\phi(x),\\sigma^2(x))\\end{align}$$\n",
    "\n",
    "w为后验均值m,方差为$\\sigma^2(x)=(\\beta^*)^{-1}+\\phi(x)^T\\sigma\\phi(x)$\n",
    "\n",
    "RVM中先关向量的数量⽐SVM中使⽤的⽀持向量的数量少得多。对于⼀⼤类回归任务和分类任务， RVM⽣成的模型通常⽐对应的⽀持向量机⽣成的模型简洁了⼀个数量级,与SVM相⽐，这种稀疏性的增⼤并没有减⼩泛化误差。\n",
    "虽然RVM的设计对矩阵求逆,时间复杂度为$O(M^3)$，高于SVM($O(N^2)$)，但是控制模型复杂度的参数以及噪声方差自动由一次训练过程确定，不需要交叉验证"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append(r\"../\")\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from prml.kernel import (\n",
    "    RBF,\n",
    "    PolynomialKernel,\n",
    "    SupportVectorClassifier,\n",
    "    RelevanceVectorRegressor,\n",
    "    RelevanceVectorClassifier\n",
    ")\n",
    "\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl4k8X2xz+TpHvLjoAii4rslE0EAREFBeQiILuKIgqC\netF7UUBRcPu5gIAICKi4IFw2d0VFFFRWBQRFBQGprLJTurdJ5vfHJKVAm6ZJ2qbN+TxPbJu8mXfS\n4jkz55z5HqW1RhAEQQg9LMU9AUEQBKF4EAcgCIIQoogDEARBCFHEAQiCIIQo4gAEQRBCFHEAgiAI\nIYo4AEEQhBBFHIAgCEKIIg5AEAQhRLEV9wQ8UalSJV2rVq3inoYgCEKJYfPmzce11pW9uTaoHUCt\nWrXYtGlTcU9DEAShxKCU+tvbayUEJAiCEKKIAxAEQQhRxAEIgiCEKEGdAxAEwX+ysrI4cOAA6enp\nxT0VIYBERkZSvXp1wsLCfB5DHIAglHIOHDhAXFwctWrVQilV3NMRAoDWmhMnTnDgwAFq167t8zgS\nAgohtNY4tbO4pyEUMenp6VSsWFGMfylCKUXFihX93tWJAwgBUjJTGPXFKMq/WJ7wZ8LpuqArvx75\ntbinJRQhYvxLH4H4m4oDCAEGvj+QY6nH2D5yO0njkvjXlf+i8/zOHEo6VNxTEwShGBEHUMr57ehv\nbDm8hXd6vkP1MtWJCoti5FUjubX+rczdPLe4pycIBSY2NhaAQ4cO0adPH4/XTps2jdTU1KKYVolE\nHEApZ+eJnbS8uCVh1nMrBa659Bp2nthZTLMShHNxOBwFfs/FF1/MsmXLPF4jDsAz4gBKOQ0qN2Dj\nwY1kOjLPef77v7+nYeWGxTQrIVRISEigXr163HbbbdSvX58+ffpkG+RatWoxZswYmjdvztKlS9mz\nZw9dunShRYsWtG/fnh07dgCwd+9e2rRpQ+PGjRk/fvw5Yzdq1AgwDmT06NE0atSIJk2a8OqrrzJ9\n+nQOHTpEx44d6dix4wVzq1WrFuPGjaNp06a0bNmSLVu2cNNNN3H55Zcze/bs7OsmTZrEVVddRZMm\nTZgwYUL28z179qRFixY0bNiQuXPP7qZjY2N5/PHHiY+Pp3Xr1hw5ciSwv9QAImWgpZx6lerRrkY7\n+i/rz4udXqRydGXe/PlNPv3zU34e/nNxT08oah56CLZuDeyYTZvCtGl5vrxz507efPNN2rZty913\n382sWbMYPXo0ABUrVmTLli0A3HDDDcyePZs6deqwceNGRo4cybfffsuoUaMYMWIEgwcPZubMmbne\nY+7cuSQkJLB161ZsNhsnT56kQoUKTJkyhVWrVlGpUqVc31ejRg22bt3Kww8/zF133cXatWtJT0+n\nUaNG3HfffaxYsYJdu3bx448/orWmR48efP/991x77bXMmzePChUqkJaWxlVXXcWtt95KxYoVSUlJ\noXXr1jz33HM8+uijvP766+c4rmBCdgAhwPxe86lXsR4d3u5A9anVWX9gPavuXEWV2CrFPTUhBLj0\n0ktp27YtALfffjtr1qzJfq1///4AJCcns27dOvr27UvTpk0ZPnw4hw8fBmDt2rUMHDgQgDvuuCPX\ne6xcuZLhw4djs5k1bYUKFbyaW48ePQBo3LgxV199NXFxcVSuXJmIiAhOnz7NihUrWLFiBc2aNaN5\n8+bs2LGDXbt2ATB9+vTsVf7+/fuznw8PD6d79+4AtGjRgoSEBK9/V0WN7ABCgEhbJM93ep7nOz1f\n3FMRihsPK/XC4vxyxZw/x8TEAOB0OilXrhxb89idFFYZa0REBAAWiyX7e/fPdrsdrTXjxo1j+PDh\n57xv9erVrFy5kvXr1xMdHc11112XXZMfFhaWPV+r1Yrdbi+UuQcC2QEIglCo7Nu3j/Xr1wOwcOFC\n2rVrd8E1ZcqUoXbt2ixduhQwhxa3bdsGQNu2bVm0aBEACxYsyPUenTt3Zs6cOdnG9uTJkwDExcWR\nlJTk89xvuukm5s2bR3JyMgAHDx7k6NGjJCYmUr58eaKjo9mxYwcbNmzw+R7FiTgAQRAKlbp16zJz\n5kzq16/PqVOnGDFiRK7XLViwgDfffJP4+HgaNmzIxx9/DMArr7zCzJkzady4MQcPHsz1vffccw81\natSgSZMmxMfHs3DhQgCGDRtGly5dck0Ce8ONN97IoEGDspPQffr0ISkpiS5dumC326lfvz5jx46l\ndevWPo1f3CitdXHPIU9atmyppSGMIPjHH3/8Qf369Yvl3gkJCXTv3p3t27cXy/1LO7n9bZVSm7XW\nLb15f0B2AEqpeUqpo0qpXP/KyjBdKbVbKfWLUqp5IO4rCIIg+E6gQkBvA108vN4VqON6DANeC9B9\nBUEIYmrVqiWr/yAmIFVAWuvvlVK1PFxyC/CuNvGmDUqpckqpalrrw4G4v1CIOJ3g1K6vTtBATkVR\ni8X1UGe/FwShRFBUZaCXAPtz/HzA9Zw4gGBDa7DbIcsOmXZwOEAp87y7FE+5/nNO/kgb52C1QHgY\nhIWBzXr2PYIgBB1Bdw5AKTUMEyaiRo0axTybEEFrsDsgI8MYfTBG3uIy5gXBqSE9E9IyjPGPDIfw\ncOMMBEEIKorKARwELs3xc3XXcxegtZ4LzAVTBVT4UwthnE7IzDLG2uk0Bt/fVbtFgcX1z0rncAZh\nNoiKNF8FQQgKiipg+wkw2FUN1BpIlPh/MeJ0Qlo6nE6C1HRjtMMLIWSjlDH44WHmnmeSISnF7DaE\nkOH06dPMmjWrwO/r1q0bp0+f9njNk08+ycqVK32dWsgTkOWYUup/wHVAJaXUAWACEAagtZ4NLAe6\nAbuBVGBIIO4rFJDsFXm6ideHFWGM3mo1D7sDEpMgMgKiIiRpHGT8cewPXtn4Cn8c/4P6lerzUOuH\nqFepnl9juh3AyJEjz3nebrdna/fkxvLly/Md++mnn/ZrbqFOQP7v01oP1FpX01qHaa2ra63f1FrP\ndhl/tOF+rfXlWuvGWms53VWUaA1ZWWbFn5ZuVvrhtuJJ0NqsZleQkQmJySbZLAQFGw9s5Nq3r+WS\nuEuY2GEil8RdQvu32rPxwEa/xh07dix79uyhadOmXHXVVbRv354ePXrQoEEDIG9Z5Vq1anH8+HES\nEhKoX78+9957Lw0bNuTGG28kLS0NgLvuuiu7J0CtWrWYMGECzZs3p3Hjxtly0seOHaNz5840bNiQ\ne+65h5o1a3L8+HG/PlNpQZZfpR2HE5JT4UyKCfWEFZPhz4k7NGRRJiyUmnZeRZFQHDz27WNM7jyZ\nJzo8QcfaHXmiwxNM6jyJx7993K9xX3jhBS6//HK2bt3KpEmT2LJlC6+88gp//vknAPPmzWPz5s1s\n2rSJ6dOnc+LEiQvG2LVrF/fffz+//fYb5cqV4/3338/1XpUqVWLLli2MGDGCyZMnA/DUU09x/fXX\n89tvv9GnTx/27dvn1+cpTYgDKK24wz2JSaasMzws+MItFotxBGkZJjfgcOb/HqHQ+C7hO/o36n/O\nc/0b9md1wuqA3qdVq1bUrl07++e8ZJVzUrt2bZo2bQp4llju3bv3BdesWbOGAQMGANClSxfKly8f\nwE9TsgkyiyAEBKdr1Z+SauLuHuKsxY5yJaAdzrPOSigWLoq5iL2n9p7z3N7TewPeN8ItAQ3nyipv\n27aNZs2aZcsq5ySnVLMniWX3dcEuwxwsiAMobWS6Yv1Z7lV/CTmIZbOaQ2SJySY/IBQ597W8j39/\n+W9Op5vKm9Pppxn15SiGtxiezzs940mSuShkldu2bcuSJUsAWLFiBadOnQr4PUoqQbw0FAqE1qak\nMz3DGNMgC/ccOgQffQR79pipNW8O3btDXFyOiywWsCmze3FqUyUkFBmPtX+MUV+MovYrtbmy4pX8\neeJPBjUaxGPtH/Nr3IoVK9K2bVsaNWpEVFQUVaqc3VF06dKF2bNnU79+ferWrVsossoTJkxg4MCB\nzJ8/nzZt2lC1alXizvmHF7qIHHRpwB3yybIHR5I3B04nzJhhjP/NNxvDn5UFq1fDuvUwZgx0PV9G\nUGvzWaIjTbloEH2ekkhB5aCPphxlz8k9XF7hci6KuagQZ1Y0ZGRkYLVasdlsrF+/nhEjRuTZeayk\n4a8ctOwASjpZdmP80QWXbSgCZs+GzVtg8bIs1h5dzmf7fiDSGkG3+7ox5O5rePBBRUw0XHttjje5\nq4RS040ziIoUJ1CEXBRzUakw/G727dtHv379cDqdhIeH8/rrrxf3lIIGcQAlmfQMSEkHmwUswae1\nc+oULF4Cy953MHHDw9iddnrV60VqVgpT10+h02W/MGHCCKZNhfbtz7PxbieQlmEOrUWLExB8o06d\nOvz888/FPY2gRBxASSRnvD/IQj45+Xw5dOwI2xJXk5yZzLwe87C6chMda1/PrYt707N3L+yOqvz6\nKzRpct4AbieQnmGS2VGRRf8hBKEUE1yZQiF/3PH+IDf+APv+hoYN4KeDP9Llii7Zxh+gfGQ5Wldv\nzebDm2lQH/btz2OQnOGgtIyimbgghAiyAyhJOJyQnAKOIor3Hz4Mf+6Eg4fg4AE4c8b0B3A4IDIS\nLqoCVS6CWrWhUSPzXA7CIyA1FcpElOV4yrELhj+acoyykWVIS89HJDTbCaSdlZgWBMFvxAGUFOwO\nc1oWjIhbYbB3L3yxHNb8AFu2mNrNnERFnRV1S0015TxurFaoXx/aXAPdbob27WndOpzXXoNne9zM\nPZ8O5cYrbqJuxSvRWvPVnhUcSjpIgzKteWITjBuXz9zcTiAl1ZwXEFnpgGO3w8cfw8KFcPw4VK4M\ngwZBjx7BfZZQ8B35s5YEsuxGy8eqjKENJAcOwHvz4cMP4I8/zHO1L4O27aBFC2jQAKpfCpdcYhyA\nG61NlvfIP7BrF2z9Gbb8bMZ6fS7ExXFN9x58fPAe/trckrHtxjDy85HUKFuDlMxkMp1ZTO0yjXff\nDqNtW6hYwYu5KmXOOJxJgbKx0mQmgOzebcp0q1SBoUOhZk34+2+YNs045+XL4fLLC3cO1113HZMn\nT6ZlS68qGEs0CQkJrFu3jkGDBhXrPMQBBDsZmZCc5qr0CVDKRmv46itjqFd+bZ5r1w6efxG6doXL\nLst/DKWgQgXzqN8Aetxink9LM0X+n3+G5YP3mZS8gD9vb8rxAQ+ybMIn7E77jUhbJBfRiHfeVGzc\nCG+8UYC5Wyxg1WY3VDY26A68lUROnoTOneHRR2HEiHNfu/NOmDULOnUym0J/ZHS01mitscjfjISE\nBBYuXFggB5CffLYvyF8imEnPMAnfQJ3stdth6RJo2wb69YHftsOjY+CX7fD5F3D//d4Zf09ERRkn\nMmMm7NwFU6ZSq2o618wbyum6Hfjp3n+YOqYR/fsrHA54c54PRsVqBTQkpYqKaAB44w1Thnu+8Xcz\ncqRZHxTIUbtISEigbt26DB48mEaNGrF//35WrFhBmzZtaN68OX379iU5OfmC9+V2zZdffknfvn2z\nr1m9ejXdu3cHYMSIEbRs2ZKGDRsyYcKE7GvykohOTk5myJAhNG7cmCZNmmSri+Y3tx07dtCqVatz\nPl/jxo0B2Lx5Mx06dKBFixbcdNNNHD5sel7t3r2bTp06ER8fT/PmzdmzZw9jx47lhx9+oGnTpkyd\nOpX09PTs+TRr1oxVq1YB8Pbbb9OjRw+uv/56brjhhoL/AfLD7ZWD8dGiRQsdkjidWqemaX3slNaJ\nSVqfSfbvkZik9eKlWtetqzWYr3Ne1/rEKf/H9uZx+ozWb7+js2rX0Rr0mSbX6OTvNvk/7onTWien\nFPdfK+j5/fffPb5ep47WGzZ4HmPDBq2vvLLg9967d69WSun169drrbU+duyYbt++vU5OTtZaa/3C\nCy/op556SmutdYcOHfRPP/2U5zVZWVn60ksvzX7+vvvu0/Pnz9daa33ixAmttdZ2u1136NBBb9u2\nTWutdc2aNfX06dO11lrPnDlTDx06VGut9aOPPqpHjRqVPc+TJ096nFtO4uPj9V9//ZV9zTPPPKMz\nMzN1mzZt9NGjR7XWWi9atEgPGTJEa611q1at9AcffKC11jotLU2npKToVatW6Ztvvjl7zMmTJ2df\n/8cff+hLL71Up6Wl6bfeektfcskl2Z/vfHL72wKbtJc2VkJAwUbOGv9ANG355RcYNwZ++AEuvwLe\nnW/CNUW5DbdYoPet2HrcAvPfJW7CBOjUBv472jwifNT8sVmN5LXVJpVBfvDXX0aiwxPNmpnrfKFm\nzZrZGj8bNmzg999/p23btgBkZmbSpk2bc67P6xqbzUaXLl349NNP6dOnD59//jkvvfQSAEuWLGHu\n3LnY7XYOHz7M77//ThPXwZKcEtEffPABACtXrmTRokXZ9yxfvjyfffZZvnMD6NevH4sXL2bs2LEs\nXryYxYsXs3PnTrZv307nzp0BcDgcVKtWjaSkJA4ePEivXr0AiIzM/SzLmjVrePDBBwGoV68eNWvW\nzO6X0LlzZypU8CZJVnDEAQQTWptSx/RM/2v8k5Ph/56D12ZBuXIw+WUYcjeEFaNchM1m5tC9O4wd\nCy88D599Cm+/A1fWLfh4OSuDbBYpVfGRqChITIRKlfK+JjHx3BqAgpBT/llrTefOnfnf//6X5/We\nrhkwYAAzZsygQoUKtGzZkri4OPbu3cvkyZP56aefKF++PHfdddc5ktLeSkR7MzeA/v3707dvX3r3\n7o1Sijp16vDrr7/SsGFD1q9ff861eamgFoScv79AIzmAYEFrSAmQ8V+xAlq1hBmvmizelq0wbHjx\nGv+cVL7IBP+XLDNnDa5tDwsW+DaWclVGJaWaQ3JCgbn5ZsjH5vG//5nr/KV169asXbuW3bt3A5CS\nkpK90vXmmg4dOrBlyxZef/317CYvZ86cISYmhrJly3LkyBG++OKLfOfRuXNnZs6cmf3zqVOnvJob\nwOWXX47VauWZZ56hf3/TQKdu3bocO3Ys2wFkZWXx22+/ERcXR/Xq1fnoo48AI0yXmpp6gUR2+/bt\nWeD6f+DPP/9k37591K3rw6KogIgDCAa0NsnejExzwMtX45+cDKMehD69jc7y1yth2nT/SjcKky5d\nYO16U246Yjg8cD9k+tALwGo560AlKVxgHngAJk2CI0dyf/3IEZg82VznL5UrV+btt99m4MCBNGnS\nhDZt2mQnZr25xmq10r17d7744ovsBHB8fDzNmjWjXr16DBo0KDuE44nx48dz6tQpGjVqRHx8PKtW\nrfJqbm769+/Pe++9R79+/QAIDw9n2bJljBkzhvj4eJo2bcq6desAmD9/PtOnT6dJkyZcc801/PPP\nPzRp0gSr1Up8fDxTp05l5MiROJ1OGjduTP/+/Xn77bfPaYJTWIgcdHGjtTndm+Xw73DTjz/CPUPh\n7wQY9RA8Pt732HpR43DAc8/C5Emm3GT+AqhYseDjZGZBTJSRkBay8UYO+tln4d13YepU45etVvNn\n+fJLePhhs5F83L/WwEIhIHLQJZlsHX+HSfj6Osar0+Gpieaw1hdfwjX5r4CCCqsVnpwAdevB/SPg\nhuth6TKoU6dg44TZXOqoQd4GMwgZPx7q1oWJE2HYMKhe3ZwRvOQSeO45yFF9KZQi5P+S4sLpNHFr\nhx/G/8QJuG84fPUl3HILvDrTJHxLKv37Q62aMHAA3NQZPv4UXDXWXqGUCQclp0IZOSRWUPr2NY9d\nu8w/rUqV4IorintWQmEi/4cUB06nOcnq8CPs8+uvcN21sOpbU+Hz7nsl2/i7ubo1fLXChK9u7gob\nNxbs/VaLEctLu7CxuOAddepA69Zi/EMBcQBFjdNptGwcTt+N/7Kl0Ol6I8b25QpT4RPEstAFps6V\n8NXXUKEi9OwB339XsPeHuc4HZGblf60ghDDiAIoSt/HXPhp/pxMmToC7h5iTOd//AKVVOKtGDbMT\nqFED+veDDevzf48bt2hcspSGCoInxAEUFY4cxt+XBGVKCtxxG0x5Ge6+Gz75zOjxl2aqVDGfs9rF\n0OdWo0bmLe74v5SGCkKeiAMoChxOSEo2hsgX43/4MHS9CT77zCh2Tn0FwkNE+qBKFfjkU5Pf6NUT\nfv/N+/eG2UwYSEJB55KYDCdOB+6ReKGYW2ESGxsLwKFDh+jTp4/Ha6dNm0ZqaqrXY3/00Uf8/vvv\nub6WkJBAo0aNvJ9oCUAcQGHjcMCZZNPY3Bf9+h1/wA0dTWnGoiVGsbO44v3FtZKuXh0+/dzU9/fu\nBQcPev9em83sAiQUdBa73Rw4DNTDg7yCtzgcjgK/5+KLL2bZsmUerwmkAyiNiAMoTBwOE/YB34z/\n2jVwY2eT7P3iKyOzXNhobXYsWXbX6tl+9pH9XJbr+xzPZdnN+wrLSdSuDe9/CElJ5qRzYqJ377O4\nnKWEgoqFhIQE6tWrx2233Ub9+vXp06dPtkGuVasWY8aMoXnz5ixdupQ9e/bQpUsXWrRoQfv27bNP\n4e7du5c2bdrQuHFjxo8ff87Y7hW5w+Fg9OjRNGrUiCZNmvDqq68yffp0Dh06RMeOHenYseMFcxs7\ndiwNGjSgSZMmjB49mnXr1vHJJ5/wyCOP0LRpU/bs2cPmzZuJj48nPj7+HOmI0oKcAygs3MZf4VsX\nr48/gqF3m9ZMH3xkvhYWTm3mq7UrgWqDiHBTUmmxGCOq1Nmdh9uQam3eq52mZWWW3XzVmKWF1RrY\n3UqjRuaUcJ/eJh+y7APvQmE5Q0ERIRI6CyJ27tzJm2++Sdu2bbn77ruZNWsWo0ePBqBixYpsceV2\nbrjhBmbPnk2dOnXYuHEjI0eO5Ntvv2XUqFGMGDGCwYMH52mE586dS0JCAlu3bsVms3Hy5EkqVKjA\nlClTWLVqFZXOU7o7ceIEH374ITt27EApxenTpylXrhw9evSge/fu2aGlJk2aMGPGDK699loeeeSR\nQvwtFQ+yAygM7A4TF/XV+L81DwbfAc2aw9ffFI7x19ps3TOzTHgkMtx02CpfBsrEQFSE2d67m9Hk\nNORuZ2CxmNfDwiAq0hy+cr8/PNz8HjKzzM4gUFx/PcyYZbqOPfyQ96t6CQUVG5deemm2Ps/tt9/O\nmjVrsl9zi6klJyezbt06+vbtS9OmTRk+fHh2Q5W1a9cycOBAAO64445c77Fy5UqGDx+e3TErP/nk\nsmXLEhkZydChQ/nggw+Ijo6+4JrTp09z+vRprr32Wo/3LsnIDiDQ2F39ey0+9O/V2ujhPPM03HiT\n0e7P5R+mXzi1McxgjH5EWGBX6m6J5jAbREcaB5CeYb5aLYHpaTxoEOzZDZNeMieF78ujlVVOLAoc\nGKmIuAD/TgWPqPP+beX82S117HQ6KVeuHFu3bvVqDH+x2Wz8+OOPfPPNNyxbtowZM2bw7bffBvQe\nJQHZAQQSux0SU1x9a30w/o+NM8a//wD436LAGn+nNiEapxNiIqF8nBFOswWg6UxeKGVCLmXjzvbv\nzcwy4SZ/eXw83Nwdxo2F71Z79x6b1aiNSlVQkbJv375smeSFCxfSrl27C64pU6YMtWvXZunSpYDR\n5t+2bRsAbdu2zW7esiAP2fDOnTszZ86cbL3/kydPAlwgu+wmOTmZxMREunXrxtSpU7PvlfP6cuXK\nUa5cuewdS173LsmIAwgUWXZT7WOzmJVuDn49+itjvxnHbe/fxsTvnmLPyfNaKzkc8O8HYeYMs5qd\nMzdw2v3aZfgdDrMiLxdnqmmKWifHZoO4GBMmQhkj7E9S1mKBua/DlVfC4MGwd2/+73EfEAv1UJAt\nR04kEI98Spvr1q3LzJkzqV+/PqdOnWJEHs2HFyxYwJtvvkl8fDwNGzbk448/BuCVV15h5syZNG7c\nmIN5VIDdc8891KhRgyZNmhAfH8/ChQsBGDZsGF26dLkgCZyUlET37t1p0qQJ7dq1Y8qUKYBpODNp\n0iSaNWvGnj17eOutt7j//vtp2rQpOZWTDx06RLdu3bz7fQcxIgcdCDKzjKSz9cLm7Wv3reOp7yZy\nb4t7aVi5ET8e/JH5v8xnRtdXqV+5vlmRDrsXPnjfNGh/fHxgVuTaFerRGqKjTKgnWMTRtDa9D1LS\nzRLEH+XOv/6C6zpAzRomX5JHy71zyLSb8FeMjy2uShjeyEEXFgkJCXTv3p3t27cXy/1LO/7KQQfE\nIiiluiildiqldiulxuby+nVKqUSl1FbX48lA3DcoyMg0wm65GH+tNa/+9CoTOkykb4O+NKhcn7ua\n3smIliOYs2UupKfDHbcb4//MszD+icAYf6erjDPMZlb8UcWw4veEUmYXUi7u7GrU14XIZZeZHdO2\nbTDGyyqNMKvJSwSgfl0QSjJ+WwWllBWYCXQFGgADlVINcrn0B611U9fjaX/vGxSkZxpJZ5stVwOb\nbk9n3+m/aXPpuY2lr699PTv3bTEaN18shylTTRMXf9HarG6d2oRa4mICk3QtLKwWiI02K/Esh++5\nga5d4T//hbfeyr+3IeSQjZazAYVNrVq1ZPUfxARiWdgK2K21/ktrnQksAm4JwLjBTXqGaUYebjt7\n2Og8wq0RRNoiOZx0+Jzn9x/4nZfn7DXJy9fmwD33+j8f96rfXc7pT3exosS9GygbY84PZPm4Kh//\nhOkm9vAo+MOLk5zullcZPrSgLIEEc6hX8I1A/E0D4QAuAfbn+PmA67nzuUYp9YtS6gulVMMA3Ld4\n0C6t+eS0fJu3Wy0Wete/lf9b83+cSj8NwJFDu4nuN4B6e06bxui33eb/nLJyrPpjooIr3OMtNpvL\ncVl9CwnZbDDvbYiJgSF3QVpa/u8Jc50NCOQ5hSAkMjKSEydOiBMoRWitOXHiBJHe5Lw8UFTLxC1A\nDa11slKqG/ARkGu/P6XUMGAYQI0aNYpoel6iNaSmm9V/uHflk8NbDmfa+qn0WtyLWo4yjJ2ygTqH\n0lHvzIcefm6U3BU+4WEl1/DnxGKB2BhITTPhtXwc7AVUrQqz58KtvWD84/DyFM/XK2V2b6lpJlxW\nSqlevToHDhzg2LFjxT0VIYBERkZSvXp1v8bwuwpIKdUGmKi1vsn18zgArfXzHt6TALTUWh/3NHZQ\nVQFpbVYoaKKoAAAgAElEQVSLGT4YJiD5cAJhvW4lfPde1PwF/uv6OJwmhBETZWrtS1NDGIC0DGOY\nffhdM26sKaldvAS6elGql5FlQlCBKr0VhGKkqKuAfgLqKKVqK6XCgQHAJ+dNqKpyHeVTSrVy3fdE\nAO5dNLibt2dmmdV2QQ3SsaPE9u5PxF9/oxYt8d/42+3GIZWNNfHz0mb8wVQuxUSfPbxWECY+BfHx\nMOI+I6WdH2FWSQgLIYnfDkBrbQceAL4C/gCWaK1/U0rdp5S6z3VZH2C7UmobMB0YoEtKQNLdv9dd\nVllQjhyBm7uZevUlS6FTJ9/n4g75WK3G+PtTP18SiAw3VUJ2Z8GcQEQEzHvLlNnePzJ/w26xmGvS\nM/ybryCUMOQgmCey5Zz9aOTSvRscOmSMf/trfZ+L2/hHhpuDXaVx1Z8XmVmucltLwfIcc+fA6P/C\ntFfg7qGer3WL45WNC+7SWUHIhyI/CFYqsdvPKnr6YvwPHoRuXY0TeP8D/4y/u8QzJir0jD+YsFts\nlDnZ7CzAguWee6Hj9UZjac8ez9e6FU5T0/2bqyCUIMQB5EZGpjH+vqpX7t8PXbvAsaPw4cdwTVvf\n5+Jwae3HxZTeeL83RLjDQXbvY/UWC8x6zSR3R9yX/0Ez96nkLBGLE0IDcQA50dpUnyTnfbo3X/7+\nG7p1gdOn4KNP4OqrfZ+PW8unbJxZBYc6EeFnE8PeOoFLLoFJL8OG9aYyKD9skhAWQgdxAG60NmWH\n7tLDPE73euSvv0zz9sQz8PGn0NKrMFzuZNnNar9MrG/tJEsr7hxIQZxA//5GOvrZZ0xvZU9kJ4RD\n44SwENqIA4CzlT4ZWb7VnQPs+tMY/7Q0+Hw5NGvm+3yy7Mbol4m5QFpawDiByAjvZSOUgqnTjFLo\nAyPzryiyWc1CoJSfEBYEsS7u9o12h+/Gf8cfJuFrtxvj37ixb3PR2sSgw1za+SX9ZG9hoZTpbRAe\nZsTvvKFqVXj+RVi/3lQH5Te+cp0QFoRSTGhbmIxMOJNkKn18FU/79VeT8FUKln8BDXyUOcpZ5hkb\nHbrJXm9RytXRzHK2xWV+DBoEnW+EiRMgIcHztTaXJpGv4nSCUAIITQfgjvcnp4LV5nvd9+bN5pBX\nVBQs/xLq1vN9Pll2E9YIxTJPX3FrB4F34Rql4JXp5u/90CjPOQTl6umcIglhofQSeg7AHe9Py/A9\n2QumquSWf0G5svDFV3DFFb6N49bwj440DzH+BcNqMU3eHQ7vDHX16jBhInz7DSxZkv/YISQZLYQe\noeUAslyHuxwO3zR93KxaBT1vgYuqmJV/zZq+jeNe+cdEQZQYf5+x2VxNZbysDBp6D1zVCsY+Cic8\n6hGasVPTQ7uHsFBqCR0H4HBAUrJZ8fujobP8c+h7K9SuDV9+aVaUvnCO8Y/wfT6CITLC5E+8yQdY\nrfDqq5CYCI8/7vla9w4xTXSChNJH6DgA98LQn8qapUvg9tugUWP4/AuzA/AFpyvsExNtDJcQGKKj\nTNjGgxPQWpOalYajXn146GFYuABWr/I8rs1qzgV4m2wWhBJC6DgAf3njdbhnKLRuDZ98ChUq+DaO\n0yU6FhdtVqxC4FDKJIW1zjVk8/VfK7l1ya10frcTN713E291roK+7HJ4+GGjHOppXHfjGEkIC6UI\ncQD5oTVMngT/eRi6dIX3P4QyZXwby+k0xj822sgaCIHHajFhNfu5SeF1+9cxdf0UHm8/njV3r+HN\nW97kh6M/svy+TrBnN0yb6nlcm9WE7DJFJ0goPYgD8ITTabpLPf0U9B8A7y0wJZ++jmV3ivEvCiJc\nJ4VzhGzm//IeD7V+iBYXN0cpRc2yNXj2+ueYErUVR+/e8PJk2L3b87g2q0kIyy5AKCWIA8iLzEwT\n8pk1E0aMhDlzfW8Z6Db+cWL8i4zoyHPyAfsT99Og8rmH9C6Oq4bVYuX0hDGmicx//+PZuItOkFDK\nEAeQG0lJ0K8vLFtq2gu+8KLvyWOnW845WhQ9ixKlzG7LqcGpubJiHX46+OM5l+w+uQelLJStWQ+e\neBJWfQsfvO95XPcuQMpChVKAOIDz+ecfI+f83WqYOQv+81/f6/OdObT8xfgXPVarKx9gZ0izIcz8\naRaf/vkZJ9NOseHARsZ8/ShDm92NzWI1zWOaNjXNY5KS8h5TKSMdkiaNY4SSjziAnOzcATd0NLHg\nxUvhjsG+jyXGPziICIPwMBqXr8/kGyfz5e4v6bPkVl7dOJ2hzYfSr2E/c53VCi9PNR3cXnzB85jZ\nZaGiEySUbEKnJ7DdYYTf8orjf/+dqfEPD4el7/sn5+x0moNncbG+i8wJgcPpNCfALSr/UN6DD8CC\n92DtOqjfIO/rHI6z/RrkBLcQREhP4IIy/10j7VC1KnyzSox/acNicbWT9EIvaOJEiI3LPyFstZrx\npCxUKMGEtgNwOODJJ+D+kXBtB/j6G991fUCMfzATZnM1kcnnNG/FSsYJrFkDS5d6vlYSwkIJJ3Qd\nQGIiDOhvDgANHQrL3oeyZX0fT4x/8BMVYcJA+UlH33mX2QU+8TgkJ+d9nbssVNRChRJKaDqA3btN\nsveblfDyFJgyzT+BODH+JQN3KCg/6Wir1TSSP3wYJr3keUyb1VQEOUQnSCh5hJ4D+OILuP46OHEC\nPv4E7h3mXxIvu9pHjH+JwB0Kyk/YrVUruO12mPGq50by2e0jpSxUKHmEjgNwOOC5Z6F/X6hVC1Z/\nD+2v9XNM1wnfsmL8SxTuxjv5xe6fespIf4x5NP8dg7SPFEogoeEATp6E7jfDlJdh8J2wYqV/yV4w\nxt/phLIx/oWPhKLHfUo4v6qgi6rAuMdh5dfwxXLP40n7SKEEEhoOICwMjh2Haa/AjJkQGenfeA6H\nMf5lYsX4l1TCbN41kBk2DOrVg7FjPUtGS/tIoQQSGg4gLg42bvTvZK8bh8Poy5SJNQlAoeQS5UUo\nKCwMXngJEvbCzBmex5P2kUIJIyQcgNMJp5OtZPq7OLM7TGexsmL8SwUWi8kH5LcLuP566N7d9IU4\ndMjDeNI+UihZlGoHsH8/jB4NF10EjRpC795w//3w/fc+DOY2EmViTLxXKB1EhButpvx0fZ573lzz\n5BOer7NZIT1D2kcKJYJS6wB+/tlU8gFs3KjZvjuRJcvsdO8Or7xi0gFe5+vs9rO6L2L8Sx/RUeDE\n8z+I2rXhwX/DksWwYX3e1ymX3pC0jxRKAKXSAaSlQY8e8Oqr0Om+L+m1Ip76M+tx28f9+a3sy8x+\nI5P162G5h8KObLIc5n/oMjEm0SeUPqwWiInMv4zzv6Ph4ovh0Uc8x/mlfaRQQiiVFm3xYmjcGGpd\ns4k7P7qTFzu9yKH/HOKNf83lcPJhZv3yAg+NggUL8lmkZdnBZjGSzr42hBFKBhHhxnB7OtEbEwNP\nPwtbt8J78z2PJzpBQgmgVFq1pUvhrrtg+sbpjGs3jq51uqKUomJ0RZ667ilWJaziiiYnOXkKDhzI\nZQCtzeotzCbGP1RQCmKizfkOT6uCvn2hdRt4aqLRk8qL7PaRkhAWgpdSadlOn4Zq1WDPqT20qNbi\nnNdiw2O4JPZijqUeoWLFXJo/aW1W/pHh5rCQaL2HDjZr/oqhSsGLL8Hx4/DSi/mPl54hOkFC0BIQ\nB6CU6qKU2qmU2q2UGpvL60opNd31+i9KqeaBuG9eVK1q9N6aVmnK1399fc5rx1KPczD5ENWiL+Wf\nf6BCxRwvag2ZdqMaGR0lxj8UcSuGegrdNGtmzpS8Ngt2/Zn3dTl1giQhLAQhfjsApZQVmAl0BRoA\nA5VS57dS6grUcT2GAa/5e19PDB4Mc+fCQ60fZs7mOUxdP5WDSQfZduRXHv7yYfo37M+Pa2KpUweq\nVnG9yek0xj82Sox/KOPt2YAnn4ToaBg3zvN1NpsJJ0pCWAhCArEDaAXs1lr/pbXOBBYBt5x3zS3A\nu9qwASinlKoWgHvnSvfukJoK779+Bd8M/ob1B9bT4a0OzN40m571etK9ynCmTYM773S9weGScy4T\nY0IAQmgTHgbhNs9O4KIqMGYsrPgKvvrK83juhLDsAoQgIxAO4BJgf46fD7ieK+g1AcNqhc8+g3ff\nhTFDGjEkdgk/9N3N2EYzObr6VobcrRg2DNq05qwgWJlYad4uGJQyu0CdT0J4+H1wRR0YNxaPx8wl\nISwEKUGXBFZKDVNKbVJKbTp27JjP41x6KWzaZE7/PvOMUX5+8klISYXXX4devTAhH4sSUTfhQqzW\n/PsGhIfD8y/A7l0wd47n8aRxjBCEBMLqHQQuzfFzdddzBb0GAK31XGAuQMuWLf3aM0dHm26PQ4cC\nduAMEMbZZG94GMRESZmnkDtRkZCRZfJDef0buekm6HwjvPA89OtnQkO5oRQo1wnh2BjJMQlBQSAs\n309AHaVUbaVUODAA+OS8az4BBruqgVoDiVrrwwG4d8FxOk2ZZ1SEKfMU4y/khVLeJYRfeMEknZ5+\n2vN1NqtZeEhCWAgS/LZ+Wms78ADwFfAHsERr/ZtS6j6l1H2uy5YDfwG7gdeBkf7e1yec2nTwio2W\nSh/BO8LDXGJxHpxAnSthxEiY/64RofKEnBAWggilg7gyoWXLlnrTpk2BGczhMB2bYqJE0E0oGA4H\nnE4yJ8PzWjQkJkLzpnDZ5bDia8+Li0zXQcOYqMKZrxDSKKU2a61benNt6MQ/rFZR8xR8w2o1IUNP\nJ4TLloWJT8HGDbBkiefxwtyS0dJDWCheQscBCII/RHpxQvi226FZc3hyPCQn532dUkaBNFkko4Xi\nRRyAIHiDNyeELRZ4aRIcPgwvT/Y8ntUqPYSFYkccgCB4S3iYyQN4quW/+moYOAhenQ579ngeL8xm\nykLlbIBQTIgDEARvcZ8Qzk8y+qmnISLCnBDObzwRixOKEXEAglAQbF6cEK5a1egEffmFFzpBIhYn\nFB/iAAShoES5BAM9rdrvGwF16sDYMZCRjwaQzWpKlOVsgFDEiAMQhILiTgh76iEcHg4vToI9u2Hm\njPzHA+MEJBQkFCHiAATBFyLCXZU8HlbtnToZbfKXXoSDuUpfnSVMQkFC0SMOQBB8QSlzkje/Cp7n\nXzShncfzaRwDEgoKIZKSYPZs6NkTunSBBx6AbduKfh7iAATBV8JsLp0gD6GgmjXhP/+FDz6A1as8\nj5cdCpKqoNLMl19C7drw9dcwaJBm9KgsqlykuflmGDgQ0tOLbi7iAATBH6IjjbH2ZLAfehhq1YZH\nRntuHAMuxdBMCQWVUjZsMC1rP/kERr20iu2Ro/k+ZTSx109hy29nyMiAIUOKbj7iAATBH6xW0zfA\nU1loZCS8+BLs3GkayXtCqbOhIE/5BaFEMnEivPCCZrf1dd5YPY3mF8XTrU43th3ZRqeF7Zj91hl+\n+AF++aVo5hM6aqCCUFg4nZCYbLSCPPWXGNAfVq8ibcM6TleKpXJ0ZWyWPMQJs+wmxBQbLbLlpYR9\n+6BdGwfbfjrOsE/uYEr36Vxa9lKjDlsulv4fDqJltZakrnyE48fh1Vd9u4+ogQpCUeKNThCQ8fxz\nZDky+fGO6xjy8RB6/O9ffP7n8twvdlcFiVZQ6UBrDvyVSd8bk/k78U/iYssb45+DQY0GsSphFS1b\n5q8iEiikEa4gBIKcOkF5SI5P3r+EZre2oNvCjXSo9BC/N6/BIytGUyG6Am2qt77wDW6toDCbyJiX\nZJxOSEmjfFgWR0/aiI0uy9GUYzi1xpJjd7f/zH4qRlckKclEDYsC2QEIQiBwt4/MQycoKTOJFX+t\noO2kxXDllfDIaBrE1mLkVfezaPuivMe0WCApVaqCSiqZWaaZUJady+uFceSowpp4GZViKjHv53k4\nXCW/u0/s4aW1L3Fv83tZtAi6di2a6YkDEIRAYbOZTl+5hIJOpp2iXGRZysZVgilTIWEvTHqJKype\nwT/J/+Q9ptVqVpCpRVgbKPiP1iaRn5Riej+E2QgPh169YOYsxfMdX2DN/jX0WNSDUV89RMd3ruOx\n9o9h2X8ta9fCoEFFM00JAQlCIImKhIwsYwBybO+rxVYjNSuNvaf2UvvaDkYy+pVp/NrQSsPKDTyP\naXN1EHOfOxCCG7sDklON4z6vjejdQ+HBB+DlZ6ow4f63cJb5i5TURMZ3eYWP3i9Dr8dgwQKIiSma\nqUoVkCAEmvQMs/o7z1gv/X0pC35dyL9bPcjljrJUua4ruypbiftmHbUq1PY8ptNpwktl48yKUgg+\ntDZJ+9Q0E7rLI2+TkQnz3oQPP4SKFaFaZTvvfx1Lo3grEyZAq1b+TaMgVUDiAAQh0GhtykLhAmP9\n7d5vWfzbEo6mHGXwL9Dr1RUwYyYMvjP/ce12Y1TiYqQ0NNhwJXrJzLpg1Z8XWVmw809wZtipUCuW\n6jUDk+gXByAIxU2WHc4kew7ZaA3dusBvv8OmTXBRlfzHzcwyctTRUYGbq+AfWXYT8tHaGP+C4joH\nEKhKLzkHIAjFjTc6QUrBK9MhNQUefcT7cdMy5HxAMKC1CfeccR0C9MX4FzPiAAShsPBGJ+jKuvDo\nGCMWt/zz/Md0S0UkSy/hYsXhMBU+aa7kvKcT4J7QGvb9Hdi5FQBxAIJQWHijEwRGLK5hQ3j4IUhM\nzH9ciwWsypwPEOnooicj0+R4HE6zy/M1H3PkCPTtA52vh+PHAztHLxEHIAiFSWSEMRCeDHV4OMyY\nZQzChCe8G9dqBaeWLmJFidaQkmocr9VidmK+8ukncPVV8P13MPpRqFAhcPMsAOIABKEwcZ8Qzm8X\n0KIFjLwf5s2D71Z7N3aY1SSF0+SQWKFjd5hVf0YWhPsR8klOhgfuh9sGQY0a8MMauHe47+P5iTgA\nQShs3DpB+TmB8U/A5VcYA5Gc7N3YYTZIzTBnD4TAo7X53Z5JArTXJZ65snkztG8H89+F0aNh5bdQ\nt15Ap1tQxAEIQmGjlCnbdOaTEI6OhtdeM7rBT4z3fuxwmwkFeWpSLxQcp9OUd6akGZkPX8s0nU6Y\nOgU63wAZ6bD8C3hyogn9FTPiAAShKLBZISocsvLZBbRuY0JBb77hfSjIXRmUlJL/LkPwjiy7Cflk\n2f1L9B49Ar16woQnoXt3WLse2rYL7Fz9QByAIBQVkRGmXjy/yp0nnjShoPtHwpkz3o1tsZixk1Kk\nPNQfAlnb/+23cE0bWL/OnPd4Zz6ULx+4uQYAcQCCUFR42TiG6GiYPQcOHPD+gBicDVEkpUh5qC84\nHMbwp2f6V9tvt8PTT0GvW0x1z+rvYMjdQSnfIQ5AEIoSbxPCV18N/x0NCxfAxx95P77Nalax4gS8\nx53oTUw6K+fgq7E+dAi6d4PJk+D2O2D199CgYWDnG0DEAQhCUeJtQhhg7Dho1hz+/W84fNj7e9hs\n5pCSOIH8cbgTven+JXrBhHzaXQPbtsHcN2DmLLObC2LEAQhCUeNtQjgsDN54A9LTYMR9BTPmYeIE\nPOKWbk5MMiGbcD9W/Q4H/N9zJuRTuTJ89z0MGBDY+RYS4gAEoTjwNiFc50p47v/g229g5oyC3UOc\nQO64yzuTU82K3+ZHovfYUWP4X3geBgyEb1cbfacSgl8OQClVQSn1tVJql+trrilupVSCUupXpdRW\npZToOwuCtwlhgKH3wL/+ZUoJCyqP7nYCZ6Q6KHvV7+rRS3iYccK+smG9Odi1fj28OsMk7ouqlVeA\n8HcHMBb4RmtdB/jG9XNedNRaN/VWp1oQSj3hYflLRoMJTcyYBdWqwd13eScYl5MwmzF+Z0L4nIAj\n56rf4l95p9YwYwZ06woRkfDNKrjzrqCs8skPfx3ALcA7ru/fAXr6OZ4ghA5unaD8JKPB1I/Pewv2\n74dRDxZcAM5mBYUpcwylE8PuCp/T7lh/mH+6O2fOwJ13wGNj4aYuJt7fpEng5lvE+OsAqmit3eUJ\n/wB5tTTSwEql1Gal1DBPAyqlhimlNimlNh07dszP6QlCkGO1QlSUd0b56tbmkNgHH8Cc2b7dy2px\n1bpnlH4VUburK1tKmhHO8yfWD/DH73BdB/j0U3jmWVj4PyhXzr8xHe7cTPHsHvL9jSilVgJVc3np\n8Zw/aK21Uiqvf1HttNYHlVIXAV8rpXZorb/P7UKt9VxgLpiWkPnNTxBKPJHhJjbtcORfhvjQw7Bx\nIzw2Dpo1M06hIFgsEKaMUXQ4TElqCQxdeMTpNAqp6ZnG4Xlqy+ktixebnVdsLHz6GbRr7994Wptw\nnEVB2ZgLekcXFfk6AK11p7xeU0odUUpV01ofVkpVA47mMcZB19ejSqkPgVZArg5AEEIOpSA2ymjP\nWCyeDbLFAnPmQodrYfAdRk7Ym17C598vzGakjbMcEBvtn7Z9sOBO8qa65LH9OdDlJiPDONvX58I1\n18Db70LV3NbDBcCpze4kItyEAItJChr8DwF9Atzp+v5O4OPzL1BKxSil4tzfAzcC2/28ryCULmw2\nUxqa39kAMGGH9xbA6dNw152QlVXw+7mdANpIHZfkkFDO6p6UNOPMAmH89+2DLjcZ4//vUfDp5/4b\nf4cDHHbjdGOji9X4g/8O4AWgs1JqF9DJ9TNKqYuVUstd11QB1iiltgE/Ap9rrb/0876CUPqI8vJs\nAEDjxjB9BqxZA4/813fj7a6DT0kreUJyWhvnl5hsqnssyj/lzpx8/bUp8dz1p3G2zz5nDub5NVc7\noKBsnFn9BwF+ZUW01ieAG3J5/hDQzfX9X0C8P/cRhJDAYoGYKGOIw1T+hqx/f5OYnPIy1KsP943w\n7b7KZTjtDnMyNirK5CWCNTegXSGU1HQzZ6s1MHF+OHuqd/Ik06f53ffgiiv8n2+W3fxOgyzn4mda\nXBCEgBIeZlaHmVne1ao/OQF27oSxY8yp4RsuWI95j80K2uJKoGaYhvYRAVpRBwKtXS0wM84mzANl\n+MH0ZB46BL7/Hm67HV6e4r+Wj8NhdnSx0UGz6s+JSEEIQrARHWm+ehMKsljg9TegQQNTn77dz/Sa\nOzdgUaYBemKSqaYpzvyAw2E0+k+dMaEehTH8gayc+eF7aN8WfvwRZr0Gr8323/hnOQAFZYIn5HM+\n4gAEIdhwh4K8PbUbGwtLlpmvt/Yyh8UCMQd3PD0l1RjftPSiyxE4nWd78Z52OSGb1f+DXOfjcMBL\nL8K/ukNsnNHyuf0O/8Z071TCrFAmJqgrrMQBCEIwEhFujJ23p3arV4cPPoTUVOjdE06eDMw83I7A\nZjUO4LTbILvCMIHaGbjr4tPcRv+MifE7OdtDIcChqIQdP3Ki67Xw7DPQp6851duokX+DOp2QaTex\n/iCo8smP4J6dIIQyMVHmq7dKng0awv8Wwd690K8vJCcHbi5KmSqY8DATgknJ4QySU4xDyLKfjXnn\n5Ri0Nq/bHWfj+WeSzQ4jMRnS0oxuQJjL6BfCAal0ewZzXupP2RtuosyW7Uy7vQ4P9ylDWpSf+QS7\nw5zsLRtrKrqCJXfiAXEAghCs5AwFebvSbtfeaAZt3gT9+pgdQWHMK9x2Ng6f5YrRJ6WYnMGpM3Aq\n0fX1jFnNnzoDJ13PnU4yq/zkVLOrcGpXeMdmDH9hrpozM9k+vCfDn/2cctWvIOz7ddw/4yeiwqKY\n9eMs38Z0V/lYlDH+/gjNFTHiAAQhmIkIN+WDBVHx7HELzH0d1q2Dgf0hPb3w5qeU6+CVa8Xu3iXY\nXKt3q8UYdKvl7AEt93Xu7635nH4OFH/uhE7X03LpDyTf3h+16jto0JAwi40HWz3IZ7s+K/iYbuMf\nHgZlYv3rKFYMiAMQhGAnKtIYyII0denbD2a+BqtXw8ABhbMT8IRSuT+KA63hzTfMwa59+xg7pCZ6\n2ivnVPmUiyxHalYKuiA5DafTGP+YKPMoASGf8xEHIAjBjsViEooFCQUB3HYbzJhpuon17lnwPgKl\ngUOH4Nbe8PBD0OYaWL+RrJu78vGOc1VrPtrxEW0uvQblrRF3x/vLxBoJjxJo/EEOgglCySDMZhKL\naRkFO/x0x2DTperee6B7N1MpVPmiwptnsKA1LFliZDIyMuClSTBsOFgsPBj1IMM/G07C6QSaVW3G\ntiNb+e7v73jtZi8ltrPsximXiSlxIZ/zkR2AIJQUoiJNHL2gXb163wqLlsCff0LnzkbfpjRz8CAM\n6A/3DjX9edeuNzIZruRyrXI1Wdh7IdViq7HuwDouiqnCgt4Luax8bc/jZtf320qF8QdQBYp5FTEt\nW7bUmwraA1UQSjMOB5xONo6goP1sN240SWG7Hd6ZDx07Fs4ciwunE95+C558wojEjX8CRt4fGEPt\nLl2Njgz6kI9SarO3rXdlByAIJQmr1fQOyLIX/BDW1VfDqu/g4otNTmDunJIrAX0+v/wCnTvBQ6Og\naVNYvxEe/HdgjL/DYR6x0WcT8qUEcQCCUNKICIeoApaGuqlZE1ashM43wuj/wl2DTV+BksqpU/Do\nI3BtO9j7F7w2x+j2X3ZZYMbPspuDaUGs5+MP4gAEoSQSHWXq5wvoBH79FZ6cXIaeGYt5+7JncHz0\nCemt2qJ/+qmQJlpIZGXB7Negabzpjzzkbtj8s6l8CsQKPaeeT9nYoNbz8QdxAIJQElHKhCTc0gr5\noDVMmgzjHoM6dWDqKxZu+PRhvn9sBUmnnehON+AY9xikpBTB5P3A6YQP3ofWrczKP74JrFkHU6dB\n+fKBuUe2fn8ExMYEvZ6PP5TeTyYIpR2r1WupiLffMUrR77yXSXrD13nkx748uK4X2zpswPLzSjbU\nvQvrzOnQ5mrTDSvYcDrhs0+h3TWmDabVCouXwMefmu5ogcLhNNIWMdEl9nBXQRAHIAglGXdjcQ+q\noekZ8N578PTTmqfWjWbHsR08dd1TTOo8iaSMMzy44RHqrniZf1/5JXYVZiSlb/kXbN1ahB8kDzIy\nYP67cPVVMGig0Q56Y55J8nbtFlgDbXcl1svGGPmNEEAcgCCUdCIjIDzcyBDnwto1cGUdOBX+C/vO\n7LVI4wMAAArkSURBVOfFzi/R8KIG1KlwBeOvHU9cRBw/Hf+WSj3asfThDfDCi7Btm0msDr7dNEkp\nahIS4OmJ0KgB3D8SIiKM4f9pE/TrF9ga/GwxN6s52WsLnfOx4gAEoaSjFMRE5pkUPnYcataCHcd3\ncPUlV2OzWHO8VdGm+jXsOL6DmjXhSGKEqZ3f9iuMfsRoCXW6HjrdAIsWBVZi+nxOnIB33zG7j/jG\nMGUKNGsOH34MP6w1hj/Qxtlt/CPCXIe7QsskhtanFYTSisUCcTHme8e5SeHYWDh5Ai6Oq8auExee\nAt514k+qxVXj5EmjGgFA2bKm3/DvO0z2+NgxGHYPXF7blI4uWQL//OPfnB0O2LIFpk4xRv+Ky+CB\n+2FvAowdB9t/hyVLTZ/jwojFZ8f7o4KuWXtRISeBBaE0YXcYTX6bNbt65dRp6NkTPvjQwbAV/eha\npxu3N7kdq8XKZzs/ZfbmOSzsuZRBt8bx2mt5lNA7nbBxAyxdCh99CMePm+fr1IHmLeDKK6FuXah2\nMZQrB+XLmZBKVhZkZprV/eFDRqbhjz9g+68mK52UZMapX9/E9Hv2gvj4wjfGWXZzj7joUhfyKchJ\nYHEAglDayMwyzVlytFF85lnIyoQRo4/wwrrn2XJ4Cwq4smJdxrQbw9eLruC332HmDC/GdzjMyds1\nP8APPxhjfuCA9/OLi4OGjUz7xdZtoEMHqFLFp49aYNwhnzBbiWjZ6AviAAQh1EnPNM3cXU4gLc0o\nI0RHmTNTV9RPxqEdHNlXlvnzYecOmDMHKlTw8X7JybBrFxw5Yk4Wnz5ldg3h4ab5S4UKZndw8cVQ\ntWrxGF63nk9URKmTdMiJOABBEEz9Z3KaabWoFBmZ5gzVsmVw/ISJEkVGQu/e0K8/xMUW94QLEfdZ\nidjogslpl0AK4gBKV/BLEISzREYYHZvUNAizERGuGDgQBgyAxDOgnSbXWwqjIGdxh3xsNiOiVwok\nnAOJOABBKM1ERQAaUtOzw0FKQbmyxT2xIiBEQj7+IA5AEEo7UZGAyt4JhIQhdFf5lIk1n1nIFfnN\nCEIoEBUBCkgp5U7AqY2kQ3iYqe8v1fEt/xEHIAihgruTVXIa2CylzzjaXdr9MdHmZG9pdXIBRByA\nIIQSEeHG8CelmARpaUiK5lz1u/skCF4hvylBCDXCbKbJCcqjimjQ467wcbjkm2OjxfgXENkBCEIo\nYrUa8bOUdCPVUNLyAg6nMfyR4SbJXdrCWUWEOABBCFUsFlMbn2E1yWGrJfhDQu7STqtVKnwCgPz2\nBCGUUcokh202Ix2RmRWcuwGtjeFXSpK8AUQcgCAIRheiTKzREEpLA5R5rriNrNZGsllhQj0RYRLu\nCSB+/SaVUn2VUr8ppZxKqTy1J5RSXZRSO5VSu5VSY/25pyAIhYRS5rxA2ThTUZNlz7XBTJHgdJrd\niN1hWl6WizNzE+MfUPz9bW4HegPf53WBUsoKzAS6Ag2AgUqpBn7eVxCEwsJqNRU1ZePMLiAzyziD\nwhaOdId5MrPM9zHRUL6MGP5CxK8QkNb6DzBt5TzQCtittf7Lde0i4Bbgd3/uLQhCIWOzmi5jDocJ\nDaVnAtokiy2WwISHtDYVPU6nGS88zJxVCIbwUwhQFDmAS4D9OX4+AFyd18VKqWHAMIAaNWoU7swE\nQcgfq9XIKkRFmBV6eqbr/IA2RtpiAYvK32BrbR5O1wNMbD88zDxydDETioZ8HYBSaiVQNZeXHtda\nfxzoCWmt5wJzwfQDCPT4giD4iMUC4RZjrJ1Os3K3O84exnJqUBq0MoZda0CZ58wTYLFy3J7IrK2v\n83XCSqIiY7m72d30b9g/v0iCUAjk6wC01p38vMdB4NIcP1d3PScIQknF4goDhdlcktO4VveuhvTa\n/R+XM1Bmh3Ay/RSt3+3ITZffxKSuUzh45iATv5vIjuM7mHjdxGL5KKFMUYSAfgLqKKVqYwz/AGBQ\nEdxXEISiRKl8D5LN3jSbdjXaMfPmmdnPta3Rlvoz6/NgqwepGF2xsGcp5MDfMtBeSqkDQBvgc6XU\nV67nL1ZKLQfQWtuBB4CvgD+AJVrr3/ybtiAIJZH1B9bTs17Pc56rGluVplWbsuXwlmKaVejibxXQ\nh8CHuTx/COiW4+flwHJ/7iUIQsnn4tiL2Xl85znPOZwOdp/czSVlLimmWYUuknIXBKHIGN5yOFM3\nTGX9/vUApGWlMWblGOpUqEODynI8qKgRKQhBEIqM5tWaM+vmWfRf1p8waxin0k7RtkZbFvdZXNxT\nC0nEAQiCUKT0rt+bHnV7sPvkbspFlqNqbG5V5kJRIA5AEIQix2axUa9SveKeRsgjOQBBEIQQRRyA\nIAhCiCIOQBAEIUQRByAIghCiiAMQBEEIUcQBCIIghChKF3aXHz9QSh0D/g7gkJWA4wEcL5iRz1r6\nCJXPCfJZ/aGm1rqyNxcGtQMINEqpTVrrPHsXlybks5Y+QuVzgnzWokJCQIIgCCGKOABBEIQQJdQc\nwNzinkARIp+19BEqnxPksxYJIZUDEARBEM4SajsAQRAEwUWpcwBKqS5KqZ1Kqd1KqbG5vK6UUtNd\nr/+ilGpeHPMMBF581ttcn/FXpdQ6pVR8ccwzEOT3WXNcd5VSyq6U6lOU8wsk3nxWpdR1SqmtSqnf\nlFLfFfUcA4UX/4bLKqU+VUptc33WIcUxT39RSs1TSh1VSm3P4/XisUta61LzAKzAHuAyIBzYBjQ4\n75puwBeAAloDG4t73oX4Wa8Byru+71qaP2uO677FtB/tU9zzLsS/azngd6CG6+eL/r+9+wnRIY7j\nOP7+1NoiipDkT7vJvwuFcJD8OcheNjcpW3KRyFE5cHDh5iAcJLnYAxt7IhetYiMlS1vaVq1FqSVq\nHbTtx2GegzZ6ZnmeGTPzfZ12nmcOn0+z/b4z88w+m3fuJnY9BZyv/bwQ+Ay05p39L7puBzYAr/7w\nfi7rUtmuADYDQ7aHbf8AuoHOKft0Ajec6AfmSlqcddAGqNvV9mPbX2qb/cDSjDM2SprjCnAcuA18\nyjJcg6XpegDosT0CYLuofdN0NTBHkoDZJANgItuY/852H0n2P8llXSrbAFgCvPtle7T22nT3KYLp\n9jhMcoZRRHW7SloC7AMuZ5irGdIc11XAPEkPJT2X1JVZusZK0/UisBb4AAwAJ2xPZhMvU7msS/Ef\nwSpA0k6SAbAt7yxNdAE4aXsyOVkstRZgI7AbmAk8kdRv+02+sZpiD/AC2AWsAB5IemT7W76xyqFs\nA+A9sOyX7aW116a7TxGk6iFpHXAV2Gt7LKNsjZam6yagu7b4LwA6JE3YvpNNxIZJ03UUGLM9DoxL\n6gPWA0UbAGm6HgLOOblRPiTpLbAGeJpNxMzksi6V7RbQM2ClpHZJrcB+oHfKPr1AV+1T963AV9sf\nsw7aAHW7SloO9AAHC352WLer7XbbbbbbgFvA0QIu/pDud/gusE1Si6RZwBZgMOOcjZCm6wjJlQ6S\nFgGrgeFMU2Yjl3WpVFcAtickHQPukzxhcM32a0lHau9fIXlCpAMYAr6TnGEUTsqup4H5wKXamfGE\nC/gFWym7lkKarrYHJd0DXgKTwFXbv3288H+W8rieBa5LGiB5Quak7cJ9S6ikm8AOYIGkUeAMMAPy\nXZfiL4FDCKGiynYLKIQQQkoxAEIIoaJiAIQQQkXFAAghhIqKARBCCBUVAyCEECoqBkAIIVRUDIAQ\nQqion5Ws0e0lGwMCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f55b3994b38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def create_toy_data(n=10):\n",
    "    x = np.linspace(0, 1, n)\n",
    "    t = np.sin(2 * np.pi * x) + np.random.normal(scale=0.1, size=n)\n",
    "    return x, t\n",
    "\n",
    "x_train, y_train = create_toy_data(n=10)\n",
    "x = np.linspace(0, 1, 100)\n",
    "\n",
    "model = RelevanceVectorRegressor(RBF(np.array([1., 20.])))\n",
    "model.fit(x_train, y_train)\n",
    "\n",
    "y, y_std = model.predict(x)\n",
    "\n",
    "plt.scatter(x_train, y_train, facecolor=\"none\", edgecolor=\"g\", label=\"training\")\n",
    "plt.scatter(model.X.ravel(), model.t, s=100, facecolor=\"none\", edgecolor=\"b\", label=\"relevance vector\")\n",
    "plt.plot(x, y, color=\"r\", label=\"predict mean\")\n",
    "plt.fill_between(x, y - y_std, y + y_std, color=\"pink\", alpha=0.2, label=\"predict std.\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.2.3 RVM for classifying"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "回到第4章的二分类模型\n",
    "$$y(x,w)=\\sigma(w^T\\phi(x))$$\n",
    "在RVM中，模型使⽤的是ARD先验，其中每个权值参数有⼀个独⽴的精度超参数。\n",
    "w后验概率分布为:\n",
    "$$\\begin{align}lnp(w|t,\\alpha)&=ln{p(t|w)p(w|\\alpha)}-lnp(t|\\alpha)\\\\\n",
    "&=\\sum_{n=1}^N{t_nlny_n+(1-t_n)ln(1-y_n)}-\\frac{1}{2}w^TAw+const\\end{align}$$\n",
    "\n",
    "由第4章的IRLS方法有:\n",
    "$$\\nabla  lnp(w|t,\\alpha)=\\Phi^T(t-y)-AW$$\n",
    "$$\\nabla\\nabla lnp(w|t,\\alpha)=-(\\Phi^TB\\Phi+A)$$\n",
    "后验概率的⾼斯近似的众数，对应于⾼斯近似的均值\n",
    "$$ w^*=A^{-1}\\Phi^T（t-y)$$\n",
    "$$ \\Sigma=(\\Phi^TB\\Phi+A)^{-1}$$\n",
    "\n",
    "拉普拉斯近似来计算边缘似然函数:\n",
    "$$p(t|\\alpha)=\\int p(t|w)p(w|\\alpha)dw\\simeq p(t|w^*)p(w^*|\\alpha)(2\\pi)^{\\frac{M}{2}}|\\Sigma|^{\\frac{1}{2}}$$\n",
    "\n",
    "对$\\alpha_i$求导令其为0解出:\n",
    "$$\\alpha_i^{new}=\\frac{\\gamma_i}{(w_i^*)^2}$$\n",
    "这个形式和RVM回归时的形式一样，再把近似的对数似然函数写成:\n",
    "$$lnp(t|\\alpha)=-\\frac{1}{2}{Nln(2\\pi)+ln|C|+(\\hat t)……{-1}^TC^{-1}(\\hat t)}$$\n",
    "其中$\\hat t=\\Phi w*+B^{-1}(t-y)$\n",
    "\n",
    "$$ C=B+\\Phi A^{-1} \\Phi^T$$\n",
    "参数更新的方法与RVM回归也一样"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
